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1.  INTRODUCTION 

The  finite  element  method  has  been  used  successfully  in  recent  times  to  solve 
various  fluid-flow  problems.  The  apparent  generality  of  the  method,  and  its 
ease  of  handling  arbitrary  flow  boundaries,  led  to  its  adoption  as  a basis  for 
work  in  computational  aerodynamics.  The  long-term  aim  of  the  work  is  to 
replace  some  elements  of  wind-tunnel  testing  of  missile  aerodynamics  with 
computer  simulation  techniques  and  thereby  gain  from  the  configurational 
flexibility  and  time  saving  of  computational  methods.  However,  reduced  staff 
ceilings  have  necessitated  the  redeployment  of  resources  on  work  with  more 
immediate  benefits  and,  as  a consequence,  the  finite  element  work  is  being 
terminated  by  this  report. 

The  project  was  initiated  by  Fletcher,  who  in  reference  1 successfully  applied 
the  method  to  two-dimensional,  incompressible  inviscid  flow,  and  later  in 
reference  2 extended  the  flow  regime  to  include  the  subsonic  case.  The  present 
work  was  intended  as  an  extension  of  Fletcher's  two-dimensional  results  to 
axisymmetric  geometries,  starting  with  the  inviscid,  incompressible  case.  The 
formulation  of  the  problem  follows  that  recommended  in  reference  1,  and  the 
program  used  for  solution  has  subroutines  copied  or  adapted  from  those  used  for 
the  results  in  reference  1. 

It  should  be  noted  that  standard  programs  are  available  which  use  singularity 
distributions  to  solve  the  present  problem,  for  example  the  method  of  Albone 
(ref. 3)  which  uses  minimal  CPU  time  but  is  limited  to  'smooth'  profiles,  and  the 
much  more  general  method  of  Hess  and  Smith(ref .4) . From  the  results  of 

reference  1 it  is  doubtful  whether  a finite  element  solution  to  the  present 
problem  would  be  much  more  efficient  in  terms  of  computing  time  than  others 
already  in  use.  The  benefit  of  the  present  work  therefore,  is  dependent  on  its 
future  extension  to  more  complicated  flow  regimes,  such  as  in  reference  2,  where 
the  singularity  methods  cannot  be  applied. 


2.  APPLICATION  OF  F.E.M.  TO  AXISYM4ETRIC  POTENTIAL  FLOW 
2.1  Field  equations 

The  governing  equations  for  the  present  problem,  in  terms  of  the 
radial  and  longitudinal  velocity  components  u and  v (see  figure  1),  are 
the  continuity  equation 


1 a 

r dr 


(ru)  + 


dv 

dz 


0 


(1) 


i 


M 

i . i 


and  the  equation  of  conservation  of  vorticity 


du  dv  n 
dz  " dr  = 


(2) 


where  r and  z are  radial  and  longitudinal  coordinates  as  shown  in 
figure  1. 


Figure  1.  Coordinate  system 

To  complete  the  problem  the  boundary  conditions  are:  that  the  velocity  must 
assume  the  free-stream  value  U in  the  far-field,  and  that  the  axis  of 
symmetry  and  body  surface  must  form  a streamline  of  the  flow.  That  is. 


’ u 

v — »-U , as  R = (r2  + z2  ) 2 -►0O,  (3) 

u = 0 on  r = 0, 

and  u = v.  rj>  (z)  on  the  body  surface  r = ro  (z) . 

Equations  (1)  to  (3)  are  identical  to  those  solved  in  reference  1 except 
for  the  non-linear  term  in  (1).  This  term  requires  a modified  formulation 
from  that  used  in  reference  1. 

2.2  Finite  element  representation 

Following  the  conclusions  of  Fletcher (ref . 1)  a Galerkin  formulation  with 
quadratic  rectangular  elements  of  ’Serendipity'  type  was  chosen.  This 
formulation  is  a Method  of  Weighted  Residuals,  where  the  weighting  function 
at  each  grid  point  is  the  shape  function  at  that  point.  These  shape 
functions,  which  are  given  in  Appendix  I,  are  chosen  so  that  the  isoparametric 
transformation 

* 

x = ^ Ni(f  ,t})Xj  where  x = r or  z,  (4) 

i=l 

transforms  axisymmetric  elements  with  nodes  (r^,  z^)  in  the  (r*z)  plane 
into  square  elements  with  nodes  r\/)  in  the  (£  ,rj)  plane,  where 

= *1  and  r\ . = ±1.  The  transformation  relates  the  shape  functions  in 
the  two  planes  as 

N.  ({ ,V)  = M.(r,z),  (5) 

where  the  point  ({  ,t?)  transforms  into  (r,z)  undeT  (4). 
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Full  details  of  this  finite  element  representation  and  the  associated 
shape  functions  N^(£,tj)  can  be  found  in  reference  1.  Its  application  to 

the  present  problem  can  proceed  in  several  ways,  as  follows. 

2.2.1  Definition  of  new  flow  variable 

Multiplying  (1)  and  (2)  by  r and  defining  the  new  flow 
variable 

P = ru  ( 


gives 


dp  dv  _ _ 

dr  + r dz 


dz  r dr 


With  the  representation  in  terms  of  the  nodal  values. 


y = /_JMi(r,z)yi 


where  y = p or  v. 


and  where  / means  the  sum  over  all  contributing  nodes  i. 


the  residuals  of  (7)  and  (8)  are  given  by 


\ dM.  „ ^ \ dM. 

Ad/  1 r ^dz1  Vi 

l l 


Rj  = Adz  Mi  • r Adr  V 


The  Galerkin  formulation 


M.  R.  drdz  = 0 

J i 


for  i = 1,2 


where  S is  the  region  of  the  (r,z)  plane  occupied  by  the  flow 
field,  then  gives  two  algebraic  equations  at  each  node  point  j: 


E(a.  .P.  + b.  .v.)  = 0 

1 ]i  i i 
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(c  , .p . - d . ,v  . ) = 0, 

v j i l J i i 


where 


M.  3r  dr  dz. 


3M. 

M .r  -5— - dr  dz  , 
3 oz 


M.  5 — dr  dz 
1 o z 


M.r  5 — dr  dz. 
3 or 


Because  the  shape  functions  M^(r,z)  are  defined  to  be  zero  in 

elements  not  including  the  k**1  node,  the  integrals  in  (15)  are 

til 

actually  evaluated  only  over  the  elements  common  to  the  j and 
t h 

1 nodes.  This  is  indicated  hereafter  by  deleting  the  sub- 
script S from  the  integral  sign. 

By  transforming  the  integrals  to  the  iso-parametric  plane  as 
described  in  reference  1 and  applying  (4)  it  can  readily  be  shown 
that 


®ji  = L_j  Gj ikzk» 
k 


D . ..  ,r.  r. , 
3 lkl  k 1* 


c . . = 

Ji 


G. ..  r. 
31k  k 


D.  ...z.r, , 
3 lkl  k 1 


where 


G...  = 
31k 


r dN.  3N  3N.  3N. 

Nj  IT  W ' 9^  W~  p dr) 


rr  r SN.  3N  3N.  3N.-i 

Djiu-  JJVilv1  irj  * *>■ 
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2.2.2  Explicit  (u,v)  formulation 
The  representation 


M.(r,z)y. 


with  y = u or  v 


can  be  used  to  express  (2)  in  a form  identical  to  that  given  for 
the  two  dimensional  vorticity  equation  in  reference  1 and  this, 
together  with  the  substitution 


M . = r.v. 

1 IX 


in  (13),  gives  (1)  and  (2)  in  terms  of  the  unknown  nodal  velocity 
components  u^,  v^,  as 


(a . . r . u . + b . . v . ) = 0 

ji  i i ji  l 


(c . .u.  - a . .v. ) = 0, 

Ji  i Ji  i 


where  a..,  b..  and  c..  are  as  defined  in  (15)  and  (16).  Note 
ji  Ji  Ji 

that  this  formulation  is  somewhat  inconsistent,  in  that  (21) 
assumes  the  representation  (18)  while  (20)  assumes  (9).  This 
choice  is  determined  by  the  necessity  to  make  the  resulting 
integrals  amenable  to  the  iso-parametric  transformation. 

An  alternative  formulation  which  may  be  more  consistent  in 
this  respect  is  derived  by  substituting 


= 9 , W • w 


into  (7)  which  gives 


) ()  a..,r.u.  + b..v.)  = 0, 

Lj  L-j  Jlk  k 1 J1  1 


i k 


where 


■jik  * JjHihwirdz 

yff„  f „ ii. 

= LU  Nj  L Nk  L3$  3t?  3t?  3£  J 

1 

3N.  3N,  3N.  3 N.  “1  ') 

+ NiL¥~  "3»T  ' W 3*  J J ^ dT?’zl 
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0 


v 

, (Dj ilk  + °jkli)  Z1 

1 


with  defined  in  (17).  The  other  equation,  (21),  remains 

unchanged. 

2.3  Boundary  conditions 

The  far-field  limit  of  the  flow  field  is  taken  as  the  semi-circle 

R = (r2  + z2  ) 5 = Ro  . and  all  velocities  are  scaled  with  respect  to  the 
free-stream  velocity  U.  Hence  for  the  explicit  formulations  described 

in  2.2.2  above,  if  the  i**1  node  is  on  the  boundary  the  conditions  (3) 
become 


u . = 0 

1 

v^  = 1 J on  R = Rj , 


u^  = 0 on  the  axis  of  symmetry  r = 0, 


and  u . 

l 


v.r0  (z.) 

l i 


on  the  body  surface. 


(25) 


For  the  (M,v)  formulation  (Section  2.2.1)  an  additional  condition  may 
be  required  on  the  axis  of  symmetry,  since  the  condition  = 0 does  not 

imply  that  u^  = 0.  Consider  an  element  on  the  axis  as  shown  in  figure  2 

with  the  corresponding  iso-parametric  element. 


Figure  2.  Iso-parametric  transformation 
for  axial  element 

Using  the  shape  functions  given  in  Appendix  1 together  with  (4)  and  (9), 
we  can  write  down  the  finite-element  representation  of  H and  r on  the 
3,6,2  boundary  of  the  element  (i.e.  the  £ = 1 boundary  of  the  iso- 
parametric element)  as 


M . 

3 


r u 
s s 


**?(!  ♦ V)  Ms  + (1  - V2  ) Ms 


(26) 
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and 

rs  = k *7(1  + V)r3  + (1 


(27) 


Hence 


us^) 


h V P3  + (1  - »?)M6 
*s  t?  r3  + (1  - Tj)r6 


(28) 


I 


■ 


i 


The  axis  condition  ug(-l)  = 0 then  gives 

Pe  = k P 3 


(29) 


provided  r6  =£  T3/4  . 

Similarly  the  value  on  the  median  line  7,5  (£  = 0 in  the  iso- 
parametric element)  can  be  written  as 


xm(T?)  = k(.V2  - 1)(X3  + X4)  + *5(1  - 1?2)(X6  + xs)+^(l  + rj)x7  (30) 


where  x = P or  r.  Thus  the  radial  velocity  component  is  given  by 

Pm  k(V  - - t?)(P6+P»)  + 'j  P7 

Um^  rm  ~ '*07  " 1)  (t3+r4  )+!s(l  - 17)  (r6  +r*  ) + k r7 

and  the  axis  condition  u (-1)  = 0,  together  with  (29)  applied  at 
nodes  1 and  2,  then  gives 

P?  = (P 3 + P4  )/2. 


(31) 


(32) 


By  applying  (29)  and  (32)  the  number  of  unknown  nodal  p values  in 
each  axial  element  is  reduced  from  five  to  two.  The  actual 
representation  for  p in  the  iso-parametric  element  can  be  shown  by 
substitution  to  be 


P(£  ,*?) 


(1  ♦ v)2 
8 


P3(l  + *)  ♦ M4(l  - $) 


(33) 


that  is,  linear  in  £ but  still  quadratic  in  17.  In  the  (r,z)  plane  it 
can  be  deduced  from  (1)  and  (2)  that  u = 0(r)  near  r = 0,  which  implies 
a form  such  as 


P = 


r2f(z) 


(34) 
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in  the  axial  elements.  If  we  choose 


Tf, 

= r3/2. 

r» 

= u/2. 

(35) 

r7 

= (r3+  r4)/2, 

then  the  iso-parametric  transformation  for  r becomes  linear,  that  is 


r 


1 + V 
4 


r3  (1  ♦ £)  + r4  (1  - £ ) 


(36) 


Thus  we  can  see  from  (33),  (34)  and  (36)  that  ft  has  the  required  form 
in  the  (r,z)  plane  if  the  geometrical  conditions  (35)  are  satisfied  and 
in  addition. 


z 

= z 

rt) 

(37) 

This  last  condition  can 

be  shown 

to  be  true  only  if 

Z|  = 

z» 

= u , 

zs  = 

Z7 

% 

(38) 

and 

z2  = 

z6 

= Zj  . 

Thus  we  deduce  that. 

for  the 

G*,v) 

formulation,  the 

elements  adjacent 

to  the  axis  of  symmetry  should  be  constructed  to  satisfy  (35)  and  (38) . 
Since  the  choice  of  the  finite-element  grid  is  largely  arbitrary, 
satisfying  these  geometrical  constraints  presents  no  problems  except  for 
those  elements  containing  the  front  or  rear  stagnation  points,  in  which 
three  nodes  lie  on  the  body  surface  and  (38)  cannot  generally  be  satisfied. 
The  boundary  conditions  for  the  (/x,v)  formulation  are  then 

"i  ' °] 

Vj  = 1 J on  R = R0  , 

fli  s 0 on  r = 0 , (39) 

/*i  = vir0(zi)r0(z.)  on  r = r„(z)  , 

together  with  (29)  and  (32). 

At  those  nodes  where  one  or  both  of  the  velocity  components  are  known, 
only  one  or  none  of  the  Galerkin  equations  was  applied,  to  equalise  the 
total  numbers  of  equations  and  unknowns.  For  the  (M,v)  formulation 
either  (29)  or  (32)  replaced  one  of  the  Galerkin  equations  at  the 
appropriate  side  node  in  the  axial  elements. 

The  far-field  condition  on  v gives  rise  to  a non-zero  right-hand  side 
in  the  matrix  form  of  the  equations  of  motion.  This  can  then  be  solved 
for  the  nodal  unknowns,  in  the  present  case  by  using  the  same  sparse - 
matrix  inversion  routine  as  described  in  reference  1. 


2.4  Order  of  the  numerical  integration 

As  in  reference  1,  two-dimensional  Gaussian  quadrature  was  chosen  to 

evaluate  expressions  such  as  (17)  because  an  m**1  order  quadrature  can 
integrate  exactly,  in  an  analytical  sense,  polynomial  integrands  of  order 
(2m-l)  or  less.  The  polynomial  order  is  dependent  on  the  order  of  the 
composite  shape  functions  which  in  this  case  is  two  or  one  (see  Appendix 
1).  However  the  total  order  can  be  reduced  by  the  differencing  factor  in 
the  integrand.  Thus  for  example,  the  maximum  polynomial  order  of  the 
integrands  in  (17)  is  expected  to  be  five  and  seven  respectively,  but  could 
be  as  low  as  two  and  three. 

Fletcher (ref . 5)  has  found  that  improved  accuracy  may  result  from 
reducing  the  order  of  integration  by  one  from  that  required  to  give  an 
exact  result.  The  effect  of  this  'reduced'  integration  was  briefly 
examined  in  the  present  work,  although  it  should  be  realised  that  such 
integration  is  still  exact  for  many  of  the  lower-order  integrals.  It  is 
clear  that  complete  reduced  integration  would  require  a variable-order 
integration  routine  and  would  nullify  many  of  the  advantages  of  the  iso- 
parametric formulation,  and  for  this  reason  was  not  attempted  in  the 
present  work. 

2.5  Geometry  of  the  flow  field 

The  body  half-profile  chosen  was  an  ellipse  of  unit  half-thickness  and 
length  2,  to  which  was  added  a sinusoidal  fore-and-aft  asymmetry  of 
amplitude  A.  The  functional  form  was  then 

r (z)  = (1  - z2 /4)^  + A Sin  (y  z)  . (40) 


A polar  grid  geometry  was  used  with  the  same  number  of  nodes  on  the 
far-field  boundary  and  body  surface,  and  Nx  nodes  along  each  segment  of 
the  z-axis  between  the  body  and  the  far-field.  A small  offset  could  be 


added  to  the  angles  of  the  polar  grid  to  make  these  asymmetrical  about 
the  r-axis.  The  spacing  of  the  nodes  could  be  varied  on  the  z-axis  and 
body  surface,  to  change  the  grid  density  near  the  stagnation  points 


relative  to  that  near  the  body  mid-point  and  the  far-field. 


3.  RESULTS 

3.1  Standard  configuration 

A "standard"  configuration  was  chosen  with  far-field  radius  RQ  = 10 

and  A = -0.3  (equation  (40))  which,  with  N = 13,  N = 21  and  conditions 

x y 

(35)  and  (38)  satisfied,  gave  the  nodal  geometry  shown  in  figure  3. 
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Figure  3.  Nodal  geometry  for  13  x 21  grid 


Initially  = 17  and  = 21  were  chosen  to  take  advantage  of  the  local 

computer  priority  system  (i.e.  less  than  2 min  CPU  time  per  calculation  on 
an  IBM  370),  but  it  was  subsequently  found  for  some  cases  that  Nx  could  be 

reduced  to  13  without  significantly  changing  the  computed  surface  velocities. 
Single  precision  arithmetic  was  used,  and  at  nodes  where  one  of  the  field 
equations  was  omitted  or  replaced,  the  continuity  equation  was  retained. 

Accuracy  of  the  results  was  judged  by  comparing  surface  velocities  with 
those  from  an  iterative  method  due  to  Albonefref . 3)  in  which  the  relative 
error  is  expected  to  be  less  than  1%.  On  figures  4 and  5 the  present 
results  are  plotted  with  squares  and  triangles  and  those  from  the  Albone 
program  with  the  broken  line,  while  the  solid  line  represents  the  body 
profile. 

3.2  Results  from  the  (A*,v)  formulation 


The  (M,v)  results  shown  on  figure  4 were  calculated  using  the  standard 
configuration  described  above,  with  'exact'  integration  and  (29)  and  (32) 
applied  except  on  the  body  surface.  The  axis  elements  shown  in  figure  3 
also  satisfy  the  additional  condition 


z5  = z?  = (zi  + z2)/2  (41) 

in  the  notation  of  figure  2,  which  was  found  to  significantly  improve  the 
results. 


SURFACE  VELOCITY 


hfSRL-0056-TH 


A FORMULATION 

0 explicit  (m,u)  formulation 

METHOO  OF  ALBONE  (REF.  3) 


' 


(2)*(l-zV«)X-0  3 Sin  (J-Z) 


Figure  4.  Results  for  surface  velocity  for  optimum  configurations 

The  (P,v)  results  shown  in  figure  4 are  generally  close  to  the  Albone 
values,  however  as  a solution  they  exhibit  a number  of  undesirable 
properties.  For  one  thing  the  axial  values  of  v,  shown  for  the  upstream 
axis  in  figure  5,  do  not  follow  the  required  monotonic  decrease  from  1 at 
infinity  to  zero  at  the  stagnation  point;  the  Albone  program  does  not 
calculate  off-body  values  and  so  these  cannot  be  shown  for  comparison. 


<n 

ui 

3 0-6 

i 


A (u.,v)  FOR  MULATION 
Q EXPLICIT  (/u,v)  FORMULATION 


Figure  5.  Values  of  longitudinal  velocity  on  upstream  axis 
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Secondly,  the  results  shown  could  not  be  improved  by  changing 
computational  variables  such  as  the  values  of  and  N^,  the  choice  of 

grid  geometry,  or  the  application  of  the  various  axis  conditions. 

Rather,  any  variation  from  the  computational  configuration  described  in 
3.1  above  could  significantly  degrade  the  results  compared  to  those 
given.  However  the  (P,v)  results  of  figure  4 were  not  significantly 
affected  by  increasing  R0  to  20  provided  the  axial  grid  was  appropriately 

scaled,  by  using  double-precision  arithmetic,  or  by  reducing  the  order  of 
integration  (from  orders  4x4  and  3x3)  either  to  all  3 x 3 or  3 x 3 and 
2x2,  even  though  the  change  in  value  of  the  equation  matrix  elements 
affected  was  typically  30%. 

Other  variations  investigated  included  the  effects  of  grid  and  body 
asymmetry.  The  accuracy  of  results  for  a symmetric  elliptical  body  (A  = 0 
in  (40))  was  similar  to  that  of  the  standard  body  shape,  all  else  being 
equal.  However  a small  grid  asymmetry  seemed  necessary  for  reasonable 
results,  even  for  the  standard  asymmetrical  profile. 

3.3  Results  from  the  (u,v)  formulations 

Results  from  the  formulation  given  by  (20)  and  (21)  are  shown  by  the 

triangles  on  figure  4.  In  this  case  the  standard  configuration  was  used 

with  the  original  values  N = 17  and  N = 21,  and  it  was  found  that  the 

x y 

same  regularised  grid  as  used  for  the  (p,v)  results  was  necessary  to  avoid 
poor  results  near  the  axis  of  symmetry,  with  the  exceptions  that  (41)  was 
not  required  but  the  body  surface  nodes  in  the  axis  elements  needed  to 
satisfy 


z6  = (zj  + zj)/2 


I 1 


and  r6  = (r3  + r3)/2. 


(42) 


* A 


It  was  also  necessary  to  use  'reduced'  integration,  that  is  third  order 

for  D..,,  and  second  order  for  C...  in  (17);  neither  exact  nor  all- 
jikl  jik  v 

third-order  integration  gave  satisfactory  results.  No  explanation  for 
these  two  requirements  is  apparent.  For  the  formulation  given  by  (23) 
and  (21),  no  reasonable  results  could  be  obtained  for  any  of  the  three 
integration  routines  or  any  combination  of  the  other  computational 
variables. 

Figure  4 shows  that  the  first  explicit  formulation  gives  results  of 
similar  accuracy  to  those  of  3.2  above.  The  present  results  also  exhibit 
similar  behaviour  on  the  symmetry  axis,  as  shown  in  figure  5,  and 
generally  have  the  same  tendency  to  deteriorate  as  the  computational 
variables  are  changed,  including  in  this  case  the  order  of  integration. 
This  deterioration  is  usually  manifest  most  obviously  in  the  stagnation- 
point  elements,  in  particular  the  values  of  tangential  velocity  at  the 
two  surface  nodes  adjacent  to  the  stagnation  points.  As  already  noted, 
these  two  nodes  cannot  in  general  be  chosen  to  satisfy  (35)  and  (38) , 
unless  the  body  surface  has  infinite  slope  for  a short  distance  from  the 
axis.  The  standard  body  shape  was  modified  to  effect  this,  but  the 
resulting  equation  matrix  behaved  as  though  it  were  singular,  and  no 
solutions  could  be  obtained. 


-1 
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4. 


DISCUSSION 


The  results  described  above  indicate  that  neither  of  the  finite  element 
formulations  is  at  present  a useful  solution  technique  for  the  problem  under 
consideration.  However  Galerkin  formulations  have  been  applied  to  similar 
problems  in  two  dimensions,  for  example  references  1 and  6.  Certain  aspects 
of  the  results  already  discussed  indicate  that  the  formulation  of  the  problem 
near  the  axis  of  symmetry  may  contribute  to  the  unsatisfactory  nature  of  the 
solutions;  for  example,  the  behaviour  of  the  v component  along  the  axis  as 
shown  in  figure  5,  and  the  constraints  imposed  on  the  nodal  geometry  and  nodal 
unknowns  by  application  of  the  p-component  boundary  condition,  which  are 
apparently  also  required  in  the  (u,v)  formulation. 

Multiplying  equation  (1)  by  r to  obtain  a tractable  finite  element 
formulation,  gives  an  equation  that  is  trivially  satisfied  on  the  symmetry 
axis;  this  is  also  true  of  (2)  in  the  (p,v)  formulation.  Thus  it  might  be 
expected  that  the  Galerkin  equivalent  of  (1)  contributes  little  at  an  axial 
node  since  the  weighting  function  (i.e.  the  shape  function  for  the  particular 
node)  biases  its  effect  towards  the  axial  node.  However  in  the  (u,v) 
formulation  where  (2)  has  a non-trivial  contribution  on  the  axis,  it  was  found 
that  the  single  equation  applied  at  the  axial  nodes  had  to  be  (1)  to  give  the 
results  in  figure  4. 

In  addition  to  the  behaviour  near  r = 0,  a major  difference  between  the 
present  formulation  and  those  of  references  1 and  6 is  that  the  latter  use 
Green's  theorem  to  include  the  boundary  conditions  implicitly  in  the  field 
equations,  as  a boundary  integral  in  the  Galerkin  formulation,  while  in  the 
present  case  the  boundary  conditions  appear  only  explicitly  in  the  equation 
matrix.  Of  relevance  here  may  be  Fletcher's  comments  in  reference  5 with 
respect  to  reduced  integration,  that  results  may  be  improved  by  relaxing  the 
constraints  on  the  polynomial  form  that  the  solution  must  assume  in  each 
element.  Inclusion  of  the  boundary  conditions  in  a Galerkin  sense,  as  in 
references  1 and  6,  may  ensure  satisfactory  behaviour  without  the  constraints 
imposed  by  conditions  such  as  (29),  (32),  (35)  and  (38).  As  already  noted, 
the  'reduced'  integration  used  here  has  little  apparent  effect  on  the  (fx,v) 
results,  but  is  required  with  the  first  (u,v)  formulation  if  any  reasonable 
results  are  to  be  obtained.  Whether  a changed  integration  scheme  such  as 
complete  reduced  integration  could  improve  the  results  from  the  second  (u,v) 
formulation,  has  not  been  determined. 

No  attempt  has  been  made  to  change  the  type  or  order  of  finite  element 
used  from  "ectangular  Serendipity.  It  was  felt  unlikely  that  the  choice  of 
element  could  be  responsible  for  the  apparent  incompatibility  in  the  formula- 
tions used. 


5.  CONCLUSIONS 

A Galerkin  finite-element  method,  using  quadratic,  rectangular.  Serendipity 
elements,  has  been  applied  to  the  equations  governing  axisymmetric,  inviscid 
incompressible  flow.  Three  different  formulations,  using  the  same  representa- 
tion for  longitudinal  velocity  and  two  different  representations  for  the  radial 
velocity  component,  have  been  tried. 

The  best  results  obtained  have,  surprisingly,  shown  only  poor  agreement  with 
correct  values  and  exhibit  a number  of  undesirable  properties,  such  as  their 
dependence  on  an  appropriate  choice  of  grid  geometry.  None  of  the  formulations 
has  provided  a useful  aerodynamic  calculation  method. 
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Ro 


Ri  . Rj 
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ijk 


r,z 


r0 
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Subscripts 

i 

1,2  etc. 


NOTATION 

amplitude  of  profile  asymmetry  (equation  (40)) 

coefficients  from  integration  in  iso -parametric 
plane  (equation  (17)) 

shape  function  in  (r,z)  plane  for  it^'  node 

shape  function  in  (£,tf)  plane  for  i^  node 

number  of  nodes  on  each  segment  of  z-axis 

number  of  nodes  on  far-field  boundary  and  body 
surface 

distance  from  origin  in  (r,z)  plane,  = (r2+z2)^ 
radius  of  far-field 

equation  residuals  in  Galerkin  formulation 
region  of  (r,z)  plane  occupied  by  flow  field 
free-stream  velocity 
coefficients  in  matrix  form  of 

field  equations  (equations  (13) , (14) , (20) , (21) , (23)) 

matrix  coefficient  (equation  (23)) 

radial  and  longitudinal  coordinates  (figure  1) 

functional  form  of  body  surface  profile  (figure  1) 

radial  and  longitudinal  velocity  components  (figure  1) 

transformed  coordinates  in  iso-parametric  plane 
(figure  1) 

derived  flow  variable,  = ru 

referring  to  i**1  node 
referring  to  nodes  1,2  etc.  (figure  2) 
values  on  one  side  of  element  boundary 
values  on  median  line  of  element 
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APPENDIX  I 

ELEMENT  SHAPE  FUNCTIONS 

The  shape  function  corresponding  to  a particular  node  is  required  to  have 
the  value  one  at  that  particular  node  and  zero  at  all  other  nodes.  The 
'Serendipity'  type  form  one  family  of  the  many  functions  which  can  be  defined 
explicitly  on  the  iso-parametric  element  (see  figure  2)  to  have  these 
properties. 

For  the  quadratic  rectangular  elements  used  in  this  report,  the  shape 
functions  of  this  family  are  defined  as  follows,  where  £ . = ±1  and  jj.  = ±1: 
at  the  corner  nodes  1 1 

i.i  N.tf.n)  = »*(i  ♦ ttpd  ♦ ♦ vvi  - l); 

at  the  side  nodes  (o,i?^), 

1-2  N.(f,T?)  = *s(l  - $2)(1  ♦ r?r?i); 


at  the  side  nodes  (t^.o). 


1.3 


N.tf.t?)  = h(l  ♦ «.)(!-  V2h 


Note  that  N^  may  be  either  quadratic  or  linear  in  the  coordinates  ((•  ,rj) . 

In  the  (r,z)  plane  the  function  M^(r,z)  corresponding  to  Ni(f,rj),  and 

defined  implicitly  by  equations  (4)  and  (5)  in  Section  2.2  of  this  report, 
will  have  a different  functional  form  in  each  element  common  to  the 

it*1  node,  and  by  definition  will  be  zero  outside  those  elements. 
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